Binary black hole merger dynamics and waveforms 
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We study dynamics and radiation generation in the last few orbits and merger of a binary black 
hole system, applying recently developed techniques for simulations of moving black holes. Our 
analysis of the gravitational radiation waveforms and dynamical black hole trajectories produces a 
consistent picture for a set of simulations with black holes beginning on circular-orbit trajectories 
at a variety of initial separations. We find profound agreement at the level of 1% among the 
simulations for the last orbit, merger and ringdown. We are confident that this part of our waveform 
result accurately represents the predictions from Einstein's General Relativity for the final burst of 
gravitational radiation resulting from the merger of an astrophysical system of equal-mass non- 
spinning black holes. The simulations result in a final black hole with spin parameter a/m — 0.69. 
We also find good agreement at a level of roughly 10% for the radiation generated in the preceding 
few orbits. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw, 95.30.Sf, 97.60.Lf 



I. INTRODUCTION 

Two black holes in a binary system spiral together due 
to the emission of gravitational waves. The final merger 
stage of a binary in which the black holes have compa- 
rable masses will produce a spectacular burst of gravita- 
tional radiation and is expected to be one of the bright- 
est sources in the gravitational wave sky. Mergers of 
binaries containing two stellar black holes are important 
targets for the first-generation LIGO gravitational wave 
detectors, now operating at design sensitivity in a year- 
long science data-taking run, as well as other ground- 
based detectors such as VIRGO and GEO. Knowledge 
of the waveforms from the final merger phase is impor- 
tant to improve the detectability of these sources by such 
detectors J,, ^] . Mergers of massive black hole binaries are 
important sources for the space-based LISA, currently in 
the formulation phase. Since LISA is expected to ob- 
serve these massive black hole mergers at relatively high 
signal-to- noise ratios , comparison of the data with cal- 
culated merger waveforms should allow a test of General 
Relativity in the dynamical, nonlinear regime. 

In the early stages of the binary inspiral, the black 
holes are widely separated and the waveforms can be cal- 
culated analytically using perturbative methods. How- 
ever, waveforms from the final merger, in which the black 
holes plunge together and form a single, highly-distorted 
black hole with a common horizon, demand full 3-D nu- 
merical relativity simulations of the full Einstein equa- 
tions. This proved be a very difficult undertaking and, 
for roughly the past decade, 3-D numerical relativity sim- 
ulations have been beset by pernicious numerical insta- 
bililtics that prevented the simulation codes from running 
long enough to evolve any significant fraction of a binary 
orbit. 

Recently, however, dramatic progress has been made 
in evolving the merger of equal mass binary black holes. 
Using excision to remove the singular regions within the 



horizons and a corotating coordinate system to keep the 
black holes fixed in the numerical grid, a binary has been 
evolved through a little more than an orbit, and through 
merger 0, i^, though without being able to extract 
gravitational waves. Simulations of excised black holes 
allowed to move through the grid on one or more orbits 
and then through merger have been carried out, with 
the extraction of gravitational radiation . In addition, 
new techniques allowing the black holes to move without 
the need for excision have been developed independently 
by the authors of this paper and another research group 
[1,01; these have been appHed to study the final plunge, 
and merger of a binary, with the calculation of gravita- 
tional waves. Recently, these techniques have been ap- 
plied to evolve a binary with nonequal masses |10| and 
the last orbit and merger of an equal mass binary [Tll |. 
It is especially noteworthy that this progress is occurring 
on a broad front, by several independent groups using 
different techniques. 

Several major open questions in the area of binary 
black hole mergers center on the dependence of the re- 
sulting gravitational waveforms on the initial data. In 
order to use numerical relativity simulations to compare 
with data from gravitational wave detectors, we need to 
model astrophysically realistic binary black hole config- 
urations. For non-spinning equal mass black holes, how 
strongly do the gravitational waveforms for the merger 
depend on the initial data? What are the effects of spin 
and non-equal masses on the resulting waveforms? The 
answers to such questions can only be approached using 
an evolutionary analysis, in which different initial data 
sets are evolved and the resulting waveforms are com- 
pared. 

In this paper, we take a step towards answering these 
questions by evolving several initial data sets for non- 
spinning equal mass black holes through the final few or- 
bits, plunge, merger and ringdown. To accomplish this, 
we apply our new methods that allow puncture black 
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holes to move freely across a grid. Using adaptive mesh 
refinement, we can resolve the dynamical regions near 
the black holes (having length scales ~ M, where M is 
the total mass of the system and we use c = G = 1) 
and the outer regions where the gravitational waves are 
extracted (having length scales ~ (10— 100)7\f). Putting 
the outer boundary typically at i « 768 Af, causally 
disconnected from the dynamical regions and wavezone 
throughout most of the simulation, we can evolve stably 
for t > 800M. The initial data sets are chosen to be 
puncture data 12]. 

Our study focuses on simulations beginning from a set 
of four inspiralling black hole configurations modelled as 
described in Sec. ^1 ^^nd using techniques discussed in 
Sec. Iml In Sec. ^ we calibrate the performance of our 
numerical techniques with a resolution study of simu- 
lations from the configuration with the shortest initial 
separation. Our main results are presented in Sec. 
where we comparatively study the radiation waveforms 
and black hole trajectories. We find a close relationship 
between the trajectory information and the waveforms, 
providing a consistent picture of the black hole dynamics. 
The radiation from the final orbit, merger and ringdown 
agrees to high precision among the runs from our range 
of initial configurations, with the runs from the farthest 
initial separation producing a promising waveform ap- 
proximately through the last three orbits 



II. INITIAL DATA 

We start by setting up initial data for equal mass bi- 
nary black holes represented as "punctures" The 
metric on the initial spacelike slice is written in the form 
gij = iJ'^Sij, where i,j = 1,2,3, with conformal factor 
ip — V'BL + u. The static, singular part of the confor- 
mal factor takes the form ipBL = 1 + X]n=i "^n/2|^— rn\, 
where the n^^ black hole has mass m„ and is located at 
coordinate position r„ . The nonsingular function u is cal- 
culated by solving the Hamiltonian constraint equation 
using AMRMG T^l. 

The punctures are initially placed on the y— axis in the 
equatorial {z = 0) plane. We need to specify the indi- 
vidual puncture mass m, coordinate position Y, and mo- 
mentum P so that the black holes are on approximately 
circular orbits and have no individual spins (they are ir- 
rotational). To accomplish this, we adapt the quasicircu- 
lar initial data from the QC sequence given in Ref. 
These data are based in turn on the results of Cook (M^; 
see also ^3); '^^o uses an effective potential method that 
minimizes the binary system energy while holding the 
orbital angular momentum fixed, to produce an approx- 
imately circular orbit. 

The parameters and physical quantities for the simu- 
lations we ran are shown in Table|lJ Here, Y is the initial 
coordinate position of each puncture along the y— axis, P 
is the linear momentum of an individual puncture, and 
m is the mass of the puncture. Mq is the initial total 



ADM mass of the binary and Jo is the initial total angu- 
lar momentum. L is the initial proper separation of the 

black holes, approximated hy L = J^^ \j9vydy, where 
each limit Hi represents a point on the Schwarzschild 
horizon associated with mass at position r^. 



Run 


±Y 


±P 


m 


Mo 


Jo 


L/Mq 


Rl 


3.257 


0.133 


0.483 


0.996 


0.868 


9.9 


R2 


3.776 


0.119 


0.488 


1.001 


0.899 


11.1 


R3 


4.251 


0.109 


0.49 


1.002 


0.928 


12.1 


R4 


4.77 


0.101 


0.492 


1.003 


0.959 


13.2 



TABLE L Initial data parameters and physical quantities for 
the runs considered in this paper. 

We note that our runs Rl, R2, R3, and R4 correspond 
closely to the initial data for models QC-6, QC-7, QC- 
8, and QC-9, respectively, along the QC sequence. As 
shown in Ref. 0| , at these relatively wide initial separa- 
tions, the parameters for these initial data sets are close 
to those derived using post-Newtonian (PN) techniques; 
see in particular Figs. 26 - 28 in Ref. [lj| . 

III. SIMULATION TECHNIQUES 

We evolve the initial data with the Hcihndol codeflTj. 
which uses a conformal (BSSN) formulation of Einstein's 
evolution equations on a cell-centered numerical grid. 
The basic equations are the same as given in Ref. (17j . 
with the exception that the evolution equation for the 
BSSN variable P* has been modified for stability as sug- 
gested in ^3 . In addition, to reduce high frequency noise 
associated with refinement interfaces, we add some dissi- 
pation of the form given in . 

Mesh refinement and parallelization are implemented 
in our code with the PARAMESH package |M EH- We 
use 4*-order centered differencing for the spatial deriva- 
tives except for the advection of the shift, which is per- 
formed with 4'^-order upwinded differencing. The re- 
finement boundary interfaces are buffered with 4'^-order- 
interpolated guard cells which, at worst, may introduce 
2"'^-order errors into second derivatives. With our cur- 
rent mesh refinement implementation, our accuracy is 
limited by spatial finite differencing error from refine- 
ment interfaces. Since we do not gain much by using 
higher order time integration, we use 2"'^-order time step- 
ping via a three-step iterative Crank-Nicholson scheme. 
Even though, overall, we expect 2"'^-order convergence, 
we have found considerable advantage in using 4"^-order 
spatial differencing over 2"'^-order spatial differencing, as 
measured by the accuracy and manifest convergence of 
the Hamiltonian constraint and other quantities. 

Traditionally, puncture black holes have been evolved 
by keeping them fixed in the grid 0, . This is accom- 
plished by factoring out the singular part i/^bl and han- 
dling it analytically, while evolving only the regular parts 
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of the metric. As explained in detail in 8' , we employ in- 
stead newly developed techniques that allow the puncture 
black holes to move freely across the numerical grid. The 
singular part is not factored out; instead the entire con- 
formal factor is evolved. Initially, the binary is set up so 
that the centers of the punctures are not located at grid 
points (as in the traditional implementation). Taking 
numerical derivatives of ipBL causes an effective regular- 
ization of the puncture singularity through the inherent 
smoothing of finite differences. 

Our new approach allows the punctures to move freely 
by implementing a modified version of the Gamma- 
freezing shift vector. The Gamma-freezing condition gen- 
erally improves numerical stability by evolving the coor- 
dinates towards quiescence, in accord with the physical 
dynamics. Our modification to this gauge is tailored for 
moving punctures, as mentioned in Specifically we 
use 

dt(3^ = laB^ (1) 

and 

dtB' = dtT' - p^djT' - f]B' (2) 

for the shift /3* and 

dta = -2aK + j3'd^a (3) 

for the lapse a. The lapse equation Q is the well-known 
"1-l-log" slicing condition, modified with an advcction 
term as in Q. We use the initial gauge conditions /3 = 0, 
5 = 0, and a — ip^'^, similar to those recommended 
in We place the punctures in the z = plane and 
impose equatorial symmetry throughout. 

We use adaptive mesh refinement to produce a numer- 
ical grid having appropriate resolution in the strong-field 
dynamical regions near the black holes and in the wave 
zone. Initially, we set up the black hole binary in a nu- 
merical domain with a box-in-box refinement structure 
having an innermost refinement oi hf and subsequent 
boxes of twice that resolution. We start with two boxes 
centered on each individual black hole, and then a box 
centered on the origin that encompasses both black holes. 
Subsequent boxes centered on the origin are used to give 
up to a total of 10 refinement levels. 

During the evolution the black holes move freely across 
the grid, changing the curvature in the surrounding re- 
gion; in response, the initial grid structure is changed 
adaptively. Paramesh works on logically Cartesian, or 
structured, grids and carries out the mesh refinement on 
grid blocks. If the curvature reaches a certain threshold 
(a free parameter in our code) at one point of a block, 
that block is bisected in each coordinate direction to pro- 
duce 8 child blocks, each having half the resolution of the 
parent block. If all points in all the child blocks fall below 
the threshold, those blocks get derefined. 

The box stretching from -48M to -I-48M in the x- 
and y— directions and from to -|-48Af in the z— direction 



is fixed in place throughout all runs.^ Boxes inside this 
one can, and generally do, change adaptively as the simu- 
lation evolves; boxes outside this one maintain their orig- 
inal locations on the grid. At the initial time, there is a 
box stretching out to 24Af; as the binary evolves, this 
box generally shrinks as the black holes spiral in towards 
the center. The resolution in the region between this 
box and the next (fixed) one at ±48M is = 32hf in 
all runs. We generally extract gravitational waves on a 
sphere of radius rex = 30M. In the early stages of the 
run, this sphere intersects the next innermost box with 
resolution /iu,/2; at later times, it is completely located 
within the region having resolution 

We extract gravitational waves from our simulations 
using the Weyl tensor component 4*4. Our wave extrac- 
tion techniques are based on Misner's method and are 
2nd_ Qj.(jgp accurate. They are robust and accurate even 
when the extraction radii cross mesh refinement bound- 
aries p^ . In particular, the waveforms computed at var- 
ious extraction radii r^x are preserved up to the leading 
order 1/r scaling and show no ill effects from passing 
through one or more refinement boundaries. 



IV. CALIBRATION OF SIMULATIONS 

We have performed detailed studies of the errors 
and convergence behavior of simulations run on the 
Haindol code using fixed mesh refinement in previous 
work 0, 0, 123, 13 ■ The simulations presented here dif- 
fer from this earlier work in two important ways: they 
are carried out using an adaptive mesh structure that 
changes as the binaries evolve, and they are run for sig- 
nificantly longer durations. In this section, we discuss 
the calibration tests we have carried out to verify that 
the code produces robust, reliable results in these more 
challenging regimes. 

We performed these calibration tests on run Rl, which 
we ran at three different resolutions. The initial data 
parameters for Rl are given in TableQ]and the simulation 
parameters for these three cases are shown in Table ^ 
We use a base resolution of p = 3M/32 = 0.09375Af from 
which we reach the three resolutions used in the runs. 
Note that the medium resolution case, hf — p/3, has the 



Rl cases 


low 


medium 


high 


hf 


p/2 


p/3 


p/4 


outer boundary 


±768M 


±192M 


±768M 




186M 


332M 


291M 



For convenience, we use M = 1 to set the scale for the compu- 
tational grid and time. Note that the actual initial total ADM 
mass for each case is Mo ~ M, as given in Table |I] 
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TABLE II: Parameters for the low, medium, and high resolu- 
tion runs of model Rl, where p = 3A//32 = 0.09375M. Tsim 
is the total duration of the simulation. 

outer boundary located relatively close, at ±192M; this 
was one of the earliest runs we did. Since the simulation 
runs for a total duration Tgim = 332M, small errors from 
this outer boundary do have time to propagate in to the 
physically interesting regions of the grid. While these 
effects are generally small and have no significant impact 
on the dynamics or wave extraction, we did use a more 
distant outer boundary at ±768Af in all other runs to 
eliminate this problem. This adds only a small overhead 
to the overall cost of the simulation, due to the fixed 
mesh refinement structure used in the outer regions. 

For these three runs, we scaled the criterion for refine- 
ment and derefinement to provide as closely as possible 
the same grid structure at the same physical time in the 
simulations. Since it is generally not possible to match 
the grid structures exactly, pointwise convergence tests 
are of questionable value. Instead, we calculated the LI 
norm of the Hamiltonian constraint Cr over the grid as 
a function of time. This norm was taken over all levels 
inside the box at ±48M, which includes the wave extrac- 
tion zone. Figure n shows the LI norms for these three 
runs; note that the initial growth of the constraint vio- 
lation is brought to a halt after approximately 50M of 
evolution and then diminishes. We attribute this behav- 
ior in Ch to a gauge wave pulse that we have observed 
leaving the source region early in the simulation. The 
gauge wave has strong high frequency components and is 
thus prone to generating differencing error and reflections 
from refinement boundaries. Cr settles down somewhat 
after the gauge wave leaves the grid, as suggested by the 
plot. The curves are scaled so that, for 2'^'^— order con- 
vergence, they would he on top of each other. As Fig. 
shows, we get 2'^'^— order convergence (or slightly better) 
for the entire course of the run. 

The gravitational waves are extracted on a sphere of 
radius rex = 30Af. The top panel of Fig. [3 shows the 
I = 2, m ~ 2 component of r'i'4 for the three differ- 
ent resolution runs, where r is the estimated areal ra- 
dius of the extraction sphere [2^. The waves start out 
in phase; note in particular the agreement in the initial 
pulse around t ^ 30A/. However, as the runs proceed, 
timing differences slowly accumulate over the relatively 
long orbital time scales. By t ~ 150M, significant dif- 
ferences have accumulated in the lowest resolution case 
(solid line) compared to the other two runs; we attribute 
this to the larger inherent numerical diffusion in this low- 
est resolution run. 

A naive convergence test, in which we take differences 
directly between the waveforms at the same simulation 
times, is shown in the bottom panel of Fig.|21 these curves 
are scaled so that, for 2"*^— order convergence they would 
coincide. Comparing the top and bottom panels, it is 
clear that the differences between these curves can mainly 
be attributed to the differences in the phases of the wave- 
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FIG. 1: The LI norm of the Hamiltonian constraint violation 
is shown as a function of time for the three different resolution 
runs of Rl given in TablelH] The high resolution case is shown 
with a dashed line. The medium (bold dashes) and low (solid) 
cases are scaled so that, for 2"**— order convergence all three 
curves would lie on top of each other. The LI norm is taken 
over all levels of the grid inside 48M, including the wave ex- 
traction region. This figure indicates satisfactory convergence 
of the Hamiltonian constraint error in our simulations. 
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FIG. 2: Gravitational waveforms and a naive convergence 
test. The top panel shows the i = 2, m = 2 mode of >I'4 
for the low (solid), medium (bold dashes), and high (dashes) 
resolution runs of Rl. The lower panel shows the differences 
between these waveforms; for 2"'^— order convergence, the 
curves would lie on top of each other. Phase difi'erences be- 
tween the waveforms account for the large difi'erences shown. 
When the phases are shifted appropriately, the convergence 
of the waves is more manifest, as in Fig.|3 



forms. At several times, the phases are actually off by tt 
so that the differences between the two waves is twice as 
large as the individual waves themselves. 

A more meaningful convergence test is shown in Fig. |31 
Here, we have shifted the waveforms in time and phase 
to adjust for the timing differences using the techniques 
described below. The time axis has been relabeled so 
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FIG. 3: Time-shifted gravitational waveforms and a physi- 
cal convergence test. The labels in the top panel are as in 
Fig. |5| In the top panel, the gravitational waveforms have 
been shifted in time so that the peak amplitude of the ra- 
diation occurs at t = 0. The differences between these time- 
shifted waveforms are shown in the bottom panel; these curves 
are scaled so that they would lie on top of each other for 2"'' — 
order convergence. 

that the peak of the radiation occurs at t = 0. The top 
panel shows the / = 2, m = 2 mode of r\E'4 for the low 
(solid), medium (bold dashes), and high (dashes) resolu- 
tion cases; note that the physical parts of the waveforms 
{t > — 125M) agree beautifully. We then carried out a 
convergence test using these shifted curves; the results 
are shown in the bottom panel. Again, the curves are 
scaled so to lie on top of each other for 2"^^— order con- 
vergence. Here we see nearly perfect convergence for the 
physical parts of the waveforms; note in particular the 
much smaller vertical scale in the lower panel of Fig. |21 



hf 




= 20M 


rex = 30M 


rex = 40M 


rex = 50M 


p/2 


Erad 


0.0347 


0.0343 


0.0336 


0.0324 




J rad 


0.217 


0.218 


0.218 


0.215 


p/3 


Erad 


0.0343 


0.0345 


0.0345 


0.0343 




J rad 


0.215 


0.223 


0.226 


0.227 


p/4 


Erad 


0.0342 


0.0344 


0.0345 


0.0344 




J rad 


0.216 


0.224 


0.227 


0.228 



TABLE III: Values of the energy Erad (in units of Mo) and 
angular momentum Jrad (in units of Mq) carried away by 
gravitational radiation for the Rl runs calculated for different 
extraction radii and different resolutions. 

As we discussed in , the gravitational waveforms can 
be used to calculate the total energy Erad and total an- 
gular momentum Jrad carried away by the radiation. We 
calculate dE/dt and dJ/dt from time integrals of alH = 2 
and I = 3 waveform components using Eqs. (5.1) and 
(5.2) in Integrating dE/dt gives the total energy 

loss due to the radiation, Erad- We find that the influ- 



ence of higher modes of the waves contribute < 1% to 
the energy and even less to the angular momentum. 

The total radiated energies calculated from the wave- 
forms extracted at different radii and with different reso- 
lutions are reported in Table lTTll For lower resolution, the 
radiated energy depends somewhat on the extraction ra- 
dius, decreasing with increasing radius. For the medium 
and high resolution cases, however, the values are almost 
independent of the extraction radius. This indicates that 
the further out, lesser refined regions of the grid produce 
significant dissipation of the waves in the low resolution 
case, whereas when the resolution is high enough, the 
lesser refined regions do not have such an effect. These 
considerations indicate that the best radius to extract 
the energy is at r^x = 30M, in a region refined enough 
that the energy does not significantly dissipate in the low 
resolution run, yet far enough way from the source that 
it changes only minimally for higher resolution runs. 

As shown in Table lTTTl the radiated angular momentum 
Jrad varies by a few percent between = 30Af and 
Vex = 50M at high resolution. This observation seems 
to agree with the notion that the angular momentum 
depends more strongly on the longer wavelength parts 
of the waves which should be extracted at greater 
distances. For all the runs reported in this paper, we use 
an extraction radius of r^x — 50M to calculate Jrad- 



V. RESULTS 

In this section we comparatively analyze the results 
of simulations based on our R1-R4 black hole configura- 
tions. Recall that the initial data for each of these sim- 
ulations starts the black holes at different separations, 
with proper separations varying from 9.9Afo to 13.2Mo, 
and provided with sufficient angular momentum that the 
runs are estimated to be on initially circular orbits. 

We can think of each of these simulations as an approx- 
imate representation of the late-time portion of an ideal 
spacetime which begins with arbitrarily well-separated 
black holes on an inspiraling trajectory which asymptoti- 
cally approaches circular orbits. As such representations, 
there are several limitations which the initial data mod- 
els may have. In particular, the initial parameters will 
only approximate the ideal trajectory since the angular 
momentum as a function of radius may differ from the 
value required for an idealized circular orbit. Likewise, 
at finite radius, the ideal spiral trajectory can only be 
approximated by our initially circular configurations. 

Furthermore the manner by which we have mapped 
from these trajectory specifications to actual initial data 
values necessarily requires making some suppositions. 
For instance, our puncture data prescription, like almost 
all field prescriptions, contains no representation of prior 
radiation generated before the time at which the initial 
data are posed. Though it is generally expected that 
the significance of such limitations on the final merger 
simulations should be reduced if the black holes begin 



sufficiently far apart, there is no clear way to assess just 
how significant such effects will be on the results, includ- 
ing the gravitational waveforms, before carrying out the 
evolutions. 

Our simulations, begun with varying initial separa- 
tions, should be affected by any initial modeling error 
in varying amounts, but should agree to the degree that 
they represent the ideal astrophysical spacetime. A key 
objective in our analysis is to identify universal char- 
acteristics among the different runs which, we reason, 
are then likely to correctly represent those aspects of the 
astrophysical equal-mass non-spinning binary black hole 
merger spacetime. 

A. Overview of Simulations 

Our comparative analysis covers four simulations la- 
beled Rl to R4 in TableU We evolved them all using the 
medium resolution of hf = p/3 except for Rl, where we 
have applied the higher hf = p/A resolution. In all runs 
we used an initial grid setup and adaptive mesh refine- 
ment as described in Sec. IIIII We evolved all the runs 
to well after the wave signal had passed the extraction 
region; the actual amount of time is noted as Tsim hi 
Table Hvl For the time-slicing condition used in our sim- 
ulations, the region where the lapse satisfies the condition 
a = 0.3 corresponds roughly with the apparent horizon 
location. We thus used the moment when the two a = 0.3 
regions around the black holes merge to specify a merger 
time Tmcrgc- The number of orbits for each run, A'orbits, 
was estimated from the trajectories shown in Fig. ^ and 
is taken up to the point at which the merger occurs. 





Rl 


R2 R3 


R4 


L/Mo 


9.9 


11.1 12.1 


13.2 


hf 


p/4 


p/3 p/3 


p/3 




421M 


531M 530M 


850M 


^merger 


l&QM 


234M 396M 


513M 


Aorbits 


1.8 


2.5 3.6 


4.2 



TABLE IV: Simulation parameters and general results. 
Tmerger Is the time at which the merger occurs, starting from 
the initial time in each run. 

A graphical overview of our four simulations is pre- 
sented in Fig. ^ showing the paths traced by the black 
hole punctures on the computational domain. These were 
obtained by numerically integrating the equation of mo- 
tion Xpunc = —f3{xpunc), which analytically gives the ex- 
act trajectory of each puncture 0. The value of the shift 
at the location of the puncture f3{xpunc) was interpolated 
between grid points as required. 

For clarity, Fig.^jshows only the track of one of the two 
black holes from each simulation. We have oriented each 
trajectory according to a physical reference discussed in 
Sec. IV Bl so that they superpose at the radiation peak. 
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x/M 

FIG. 4: Paths of black holes starting from different initial 
separations. For clarity, we show only the track of one of the 
black holes from each simulation. The paths are very similar 
for approximately the last orbit, indicating that the black 
holes follow the same tracks. The point of merger (estimated 
by a single connected isosurface of a = 0.3) is indicated by as 
an asterisk in the plot. 

which occurs very near the end of the puncture trajec- 
tory. R4 has the widest initial separation and completes 
the largest number of orbits. Each of the other cases, 
after an initial transient period of approximately one or- 
bit, nearly locks on to the R4 trajectory. For the final 
orbit, all four trajectories are very nearly superposed. In 
Sec. IV CI we study the quantiative features of these tra- 
jectories in more detail. 

Fig. [SI shows one polarization component of for 
the runs Rl - R4, extracted at r^x = 30A/ and shifted in 
time and phase as described below. Notice that beyond 
about t = —50Mf the waveforms superpose sufficiently 
well that it is not possible to distinguish the curves in the 
plot. The very strong agreement among our simulations 
on this part of the wave gives us confidence that our 
waveform accurately represent the astrophysical merger- 
ringdown signal to high precision. The inset shows that 
the agreement remains generally good going back to the 
beginning of each simulation, with the R3 and R4 runs 
agreeing fairly well over some 450Af/. 

B. Gravitational Radiation 

The first step in quantitatively comparing our runs is 
to calculate the energy and angular momentum carried 
away by the gravitational radiation generated by the bi- 
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FIG. 5: Waveforms from runs Rl - R4. The figure shows 
nearly perfect agreement after t = ~50Af/. For the preceding 
500M/, shown in an inset, the waveforms agree in phase and 
amphtude within about 10% except for a brief initial pulse at 
the beginning of each run. 

nary system in our simulations. Following the discussion 
in Sec. lIVI we measure the radiation energy extracted at 
Tex = 30Af, estimating that these will be accurate within 
^ 1%. We subtract the radiation energy Erad from the 
initial mass Mq, given in Table to determine the fi- 
nal black hole mass in each simulation, Mf — Mq — Erad- 
This provides a physical scaling which we use to compare 
the R1-R4 simulations. Similarly, we measure the angu- 
lar momentum content of the radiation Jrad as extracted 
at Tex = 50M, estimating the accuracy to be within a 
few percent. Subtracting this from the initial ADM an- 
gular momentum Jq, we calculate the spin parameter of 
the final black hole a/Mf = (Jq — Jrad)/Mj. The results 
are summarized in Table fVl 





Erad/Mf J 


rad/ Mf 


a/Mf 


MQN/Mf 


aQN/Mf 


Rl 


0.0356 


0.246 


0.694 


1.005 


0.721 


R2 


0.0369 


0.272 


0.691 


1.002 


0.686 


R3 


0.0381 


0.306 


0.689 


1.004 


0.694 


RA 


0.0387 


0.325 


0.702 


1.004 


0.693 



TABLE V: Energy and angular momenta for the radiation and 
final black hole. Erad and Jrad are measured at Tgx = 30M, 
and Tex ~ 50M, respectively. Mqn and uqn are calculated 
independently from the quasi-normal fits of the ringdown 
waveforms, and agree well with the values deduced from the 
radiative losses. 

We note that the energy and angular momentum con- 
tent of the radiation is almost entirely contained in the 
^ = 2, m = ±2 spin —2- weighted spherical harmonic com- 
ponents, with other components entering at the 1% level. 
In the remainder of our waveform analysis we concentrate 
exclusively on the leading component, 

rMd,^) - r^M22) (-2^2 2(0,0) + -2>2-2(0,<^)) (4) 




-600 -400 -200 200 



FIG. 6: Amplitudes, absolute value of (complex) ^'4, of the 
waves. The curves have been shifted such that the maxima 
are all at time 0. The inset zooms into the peak showing 
the strong agreement from t = — 50M/ on. We have used 
the amplitude peak as a reference to align our simulations in 
time. 



where, for simplicity, we have suppressed retarded time 
dependence. Hereafter we suppress the multipole labels 
and refer to the leading component simply as The 
two polarization components of the radiation are repre- 
sented in the real and imaginary parts of r\E'4. We fol- 
low Refs. 10, H^, representing our waveforms by — 
Aexp{—iippoi), where the amplitude A and polarization 
phase ippoi, like r^4, are functions of time. Of course, any 
complex time-series could be represented in this form, 
but it is particularly valuable if the radiation exhibits 
a circular polarization pattern, so that A and ippoi vary 
slowly compared to the wave frequency timescale. Such 
radiation will be circularly polarized to an observer on 
the system's rotational axis, varying to linearly polariza- 
tion for an observer on the equatorial plane. We find that 
the radiation produced in our simulations shows strong 
circular polarization, which we will utilize in comparing 
the radiation from our four simulations. 

As we are interested in the degree to which our vari- 
ous runs may be taken to be different models of the same 
physical merger spacetime, we must define a physical ba- 
sis for comparing them. As expected, runs beginning 
with more separated black holes take a longer time to 
reach the point of merger. For each run, a physical ref- 
erence time is recognizable by the point at which the 
radiation reaches its peak amplitude; we define t = at 
this point for our comparisons. 

Fig. shows the wave amplitudes A{t) from all four 
runs. Through the strong radiation peak after t = 
— 50-A'f/ all four wave amplitudes show striking univer- 
sality in the compared results, with agreement among 
all runs to about 1%. This period of strong agreement 
covers roughly the last orbit, plunge and ringdown of 
the merger. Agreement within about 10% is maintained 
among the R2-R4 simulations for most of each run with 
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FIG. 7: Graviational wave phase angle vs. time. The phase is 
made to agree at time t — OMf. After t = — 50M/ the phases 
agree very well for all the runs. For the duration of the runs 
(except for a brief initial period), the waveforms from all runs 
agree in phase within about 10% of a wave cycle. 



slightly more difference in the Rl run. The smooth shape 
of the peak, lacking any sign of the several wave cycles 
spanning the peak (see Fig. [SJ is an indication of the 
circular polarization pattern discussed above. The only 
clearly non-universal feature in the wave amplitudes is a 
small and brief burst lasting about 50M/ at the begin- 
ning of each run. We interpret this burst as "spurious 
radiation" content in the initial data. As generally ex- 
pected the amplitude of the spurious radiation lessens as 
the initial separation of the black holes is increased. 

Since the black holes in our various runs (which all be- 
gin on the y-axis) undertake differing amounts of orbital 
motion before merger, we must expect differing orienta- 
tions for the systems at the point of merger. For mean- 
ingful comparsions we must rotate the data from each 
system to align them with repect to some physical ref- 
erence. As our reference, we will orient the systems so 
that the polarization phase ippoi of the radiation passes 
though zero at the moment of peak amplitude t — 0. 
We use this orientation for all figures with wave phase 
information, and in the trajectory comparison in Fig. ^ 

The polarization phase ippoi is shown in Fig. \7\ At 
time t — QAIf the phases are set to agree. However, 
they keep agreeing after that, showing that the ringdown 
frequency is the same. The plot shows generally good 
phase agreement among all runs, to within a small frac- 
tion of a wave cycle, except for a brief period at the 
beginning of each run when the radiation is dominated 
by spurious radiation associated with initial data model- 
ing error, and is not circularly polarized. For circularly 
polarized radiation it is meaningful to define an instan- 
taneous frequency for the wave, w = dippoi/dt, shown in 
Fig. |S1 which allows a more detailed comparison of the 
simulation waveforms. It can be seen that apart from the 
noise due to the smallness of the waves initially and at 
the end, all runs have the same frequency evolution from 
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FIG. 8: Waveform frequency as a function of time. The curves 
are nearly indisinguishable after t = —50Mf. 

about 60Mf before the merger. Before that the frequency 
evolution compares similarly among the simulations with 
some wavyness which may correspond to some ellipticity 
in the early part of each simulation's insprial. 

The final (constant) frequency is the frequency of the 
ringdown. Fitting an exponential decay to the ampli- 
tudes and using this final frequency we estimate, using 
the procedure in Ref. ^2^, the mass and angular mo- 
mentum of the black hole formed in the merger. These 
estimates provide a description of the final black hole as 
determined by its perturbative dynamics. The results 
are listed in Table IVl where they can be compared with 
the independent estimates obtained by subtracting the 
radiation losses from the initial values. The good agree- 
ment, to 1 — 2%, provides a measure of energy and angu- 
lar momentum conservation in each simulation. We also 
note the strong agreement among our simulations (which 
start from different separations) on the final state of the 
remnant black hole, with the measures from the R2-R4 
runs consistent with the same value, a/Mf = 0.69(±1%). 
The Rl value differs by a few percent, which may be a 
consequence of the shorter duration of this simulation, 
allowing greater sensitivity to initial transient effects. 

C. Trajectory analysis 

In this section we consider the particle-like dynamics 
of our simulations defined by the coordinate trajecto- 
ries of our black hole punctures. It is important to be 
careful in interpreting such coordinate-dependent infor- 
mation which may include non-physical gauge features. 
Nonetheless it is worth noting that the Gamma-freezing 
gauge condition applied in our simulations, to the extent 
that it approximates F* — = 0, is similar to the Dirac 
gauge applied in some post-Newtonian calculations, and 
might be expected to provide a sensible coordinate sys- 
tem at least in weak field regions. In any case, since 
we apply similar coordinate conditions in each run, the 
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FIG. 9: The coordinate separation between the punctures is 
shown for runs Rl - R4 as a function of time. Early on in 
each simulation the separation seems to drop quickly, but all 
runs track together as they approach merger. 

coordinate-based puncture trajectories provide another 
opportunity to identify universal features in our simula- 
tions that start from different initial separations. 

Refer again to Fig. ^ which shows the tracks of one 
of the punctures from each of the runs Rl - R4 as they 
spiral into the center. As the black holes descend deeper 
into the strong field, we see that their paths lock on to 
a universal trajectory that takes them into the plunge 
and subsequent merger. This behavior can also be seen 
in Fig. 1^1 which shows the coordinate separation between 
the punctures as a function of time. In this case, there is 
strong agreement among the runs after t « — 50M/. In 
the earlier part of the simulations there are clear differ- 
ences among the runs, which are suggestive of ellipticity 
in the initial orbital motion. 

One way to estimate the utility of such coordinate in- 
formation is by comparison with our much more invariant 
waveform data. We can, for instance, compare the coor- 
dinate orbital frequency with the waveform frequency 
examined above. With suitable coordinates, in the weak- 
field limit, we would expect the orbital frequency to be 
approximately equal to half the gravitational wave fre- 
quency. In Fig. IIUI we compare the coordinate orbital 
frequency with half the gravitational wave frequency 
from our R4 run. We find that, if we shift the orbital fre- 
qency data by about 33A/o, the two curves match very 
well. Despite some noise in our wave frequency, most 
features in the orbital frequency are tracked in the wave 
frequency as well. Exceptions occur at the very begin- 
ning of the simulation, where our coordinates necessarily 
start with the punctures non-moving (hence f2 = 0) be- 
fore the brief period through which the coordinate punc- 
ture velocities seem to adapt well to the physical dy- 
namics. Likewise, at late times, we see that the orbital 
puncture frequency continues to grow while the radiation 
frequency saturates at the quasinormal ringing frequency. 
This seems to correspond to the expectation that, at late 
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FIG. 10: Orbital and radiation frequencies for R4. We show 
the evolution of the orbital frequency fl, calculated from the 
coordinate motion of the punctures, compared with the wave 
frequency. The good correspondence motivates a closer exam- 
ination of the coordinate-based puncture trajectories. (Note 
that the radiation frequency is divided by two since uj — 2Q.) 



times, the coordinate motion of the punctures decouples 
from the process of radiation generation, with the punc- 
tures continuing to fall into the newly formed black hole 
while gravitational ring-down radiation is generated in 
the final black hole's pertubative potential barrier region 
somewhat outside the horizon. 

The 33Mo time shift required to realize this agreement 
can be interpreted as the time required for the gravita- 
tional radiation to propagate out to the point where it is 
extracted at r^x = 30M. We can utilize this time corre- 
spondence to roughly associate phenomena occurring in 
the strong field region of the simulations with features 
in the radiation. We have used this time-shift to com- 
pare the time at which peak radiation is generated in our 
simulations with Tmerge in Table IWI In our simulations, 
we find that the merger occurs about 12M before the 
radiation peak is generated. 

In Fig. ^] we compare the orbital frequencies of the 
puncture trajectories among our R1-R4 runs. As was the 
case with the wave frequencies we see excellent agreement 
among the runs after t = — 50M/. Here we have focused 
on the earlier region where the results are not quite as 
universal. In each run we notice that after a period of 
initial angular acceleration the orbital frequencies agree 
to within a few percent. This initial transient period is 
at least partially a coordinate effect (since f2 must begin 
at 0), but the strong similarity with the wave- frequencies 
suggests that there may be a physical basis for the de- 
screpancies as well. Particularly for R3 and R4, the shape 
of the curves suggests a slow oscillation, perhaps some el- 
lipticity in the motion. For an external comparison, we 
have included the post-Newtonian frequency evolution 
1.5PN and 2PN order, as provided in Ref. |23, positioned 
in time so that all curves agree near t = —200 A//. The 
correspondence with these PN curves is certainly strik- 
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FIG. 11: The frequencies as functions of time. Shown are 
frequencies calculated from the waves showing the agreement 
of the different runs and two curves calculated using Post- 
Newtonian approximations. 

ing, though we caution that the 2.5PN curve does not 
agree so well, with its peak frequency topping out too 
low for a useful comparison here. In subsequent work 
we intend to explore comparisons with post-Newtonian 
calculations in more detail. 



VI. DISCUSSION 

We have presented a set of numerical simulations rep- 
resenting the last few orbits and merger of an equal-mass 
non-spinning binary black hole system. Though the ini- 
tial data differ, with the black holes starting out at sep- 
arations ranging from 9M to 13Af , the calculated wave- 
forms are dominated by universal characteristics. Over 
the period beginning from about 50M before the gravi- 
tational wave peak, and covering approximately the last 
orbit, all our simulations show profound agreement, with 
differences among the waveforms at no more than about 
1%. The robustness of this part of our waveform is strong 
evidence that it accurately represents the final burst of 
gravitational radiation from an astrophysical system of 
equal-mass, non-spinning black holes, as predicted by 
Einstein's field equations. 

We also see good agreement among our simulations in 
the gravitational radiation generated in the several pre- 
ceding orbits which we have simulated. Excepting an ini- 
tial transient period about lOOAf for each run, the earlier 
portion of the waveforms show good agreement, within 
approximately 10% in phasing and amplitude. For our 
longest simulation (R4) this suggests that we have pro- 
duced a good approximation of the astrophysical wave- 
form prediction covering more than 400M before the ra- 
diation peak. We are further encouraged by good agree- 
ment over much of this period among the frequency evo- 
lution of our simulated waveforms (especially R4), with 
the second-order (2PN) post-Newtonian predictions. We 



are currently exploring the correspondence of our wave- 
forms with post-Newtonian calculations, to be developed 
in more detail in a future publication. 

That the late-time part of the simulations (the final or- 
bit and thereafter) shows stronger universality than the 
earlier part of the simulations supports idea that the late 
dynamics of the system is dominated by the strong inter- 
action of the holes, and radiative losses, which have the 
effect of reducing dependence on prior conditions. For 
the remnant black hole formed in the merger, our sim- 
ulations consistently predict a spin within about 1% of 
ajm — 0.69. 

Our simulations have employed newly developed nu- 
merical relativity techniques for evolving black holes, 
which allows the black holes to propagate accurately 
across the numerical domainj^Q. This approach does 
not require excision of the black hole interior from the 
computational domain. We have calibrated our approach 
on one of our black hole configurations, Rl, demonstrat- 
ing 2"'^-order convergence, and waveform accuracy at the 
1% level in the last orbit. For all our simulations we have 
found good energy and angular momentum conservation 
as measured by comparing the mass and spin of the fi- 
nal black hole, measured by its quasinormal ringing, with 
the expected remainder after radiative losses. In future 
work we expect to continue to refine our techniques for 
accuracy and efficiency over long-lasting simulations. 

We have also studied the trajectories traced out by the 
motion of the black holes in our numerical coordinate sys- 
tem. These trajectories suggest a coherent picture of the 
system's evolution which is qualitatively, and in some 
ways, quantitatively consistent with invariant informa- 
tion measured in the radiation waveforms. This corre- 
spondence provides some foundation for giving a tenta- 
tive physical interpretation to some coordinate-based in- 
formation in our simulations, such as the number and 
rate of inspiral orbits, and recommends further research 
toward a deeper understanding of the properties of the 
coordinate systems which we have applied in these sim- 
ulations. 

Our comparative analysis provides some insight into 
the quality of the initial data models which we have ap- 
plied. We have seen an indication of a small amount 
(with negligible energy content) of spurious radiation 
arising from initial modeling error. As generally expected 
the scale of this spurious radiation decreases for increas- 
ingly well-separated initial data. For sufficiently long- 
lasting runs this spurious radiation is well-segregated in 
time, limiting its direct significance in interpreting the 
merger radiation. Of perhaps greater concern are indi- 
cations which suggest ellipticity in the inspiral trajec- 
tories and waveforms. Two aspects of our initial data 
model may be contributing to this effect. First, the will 
be some level of mismatch between the initial separa- 
tion of the black holes and the specified initial value of 
angular momentum. We comment, in this regard, that 
other simulations we have looked at in preparation for the 
work presented here hint that the early transient part of 
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each simulation, as well as measures such as the merger 
times, may depend sensitively on small (1%) changes in 
the initial angular momentum. Secondly, as is custom- 
ary, the data we have applied have the black holes set on 
initially circular trajectories, with vanishing radial mo- 
mentum components. This may also be expected to lead 
to ellipticit y. as has been noted in binary neutron star 
simulations [23. We expect to explore these and similar 
concerns in future simulations. 

Taken together with other recent progress in numeri- 
cal relativity, these results herald a new age of numer- 
ical simulations applied to further characterize and un- 
derstand strong-field binary black hole interactions and 
merger radiation. Future applications will begin to ex- 
plore the physical parameter space of these systems to 
study, particularly, the effect on the radiation waveforms 
of individual black hole spin, and variations in the black 



hole mass ratio. 
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